Loading the packages and datasets
#Download packages
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(sp)
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(rgeoda)
## Loading required package: digest
library(spdep)
## Loading required package: spData
##
## Attaching package: 'spdep'
## The following object is masked from 'package:rgeoda':
##
## skater
#Set working directory
setwd("/Users/rebeccashi/Desktop/GGIS224/Final Project")
#load in chives
chives = st_read("chives-data-minusNU.geojson")
## Reading layer `chives-base' from data source
## `/Users/rebeccashi/Desktop/GGIS224/Final Project/chives-data-minusNU.geojson'
## using driver `GeoJSON'
## Simple feature collection with 801 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -87.94025 ymin: 41.64429 xmax: -87.52416 ymax: 42.02392
## Geodetic CRS: WGS 84
#load in variables plant and human
plant = dplyr::select(chives,"trees_n")
heat = dplyr::select(chives,"heatisl")
human = dplyr::select(chives, "acs_population")
#summary
plant_summary = summary(plant$trees_n)
head(plant_summary)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 874.000 1452.000 2135.905 2323.000 58668.000
heat_summary = summary(heat$heatisl)
head(heat_summary)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 23.00000 54.20000 71.30000 67.66205 82.90000 97.80000
human_summary = summary(human$acs_population)
head(human_summary)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2001.000 3077.000 3422.213 4500.000 19889.000
#plant_percentage
plant[is.na(plant)] <- 0
plant_sum = sum(plant$trees_n)
plant$plant_percentage = plant$trees_n/plant_sum
#human_percentage
human[is.na(human)] <- 0
human_sum = sum(human$acs_population)
human$human_percentage = human$acs_population/human_sum
#plant map
plant_map = tm_shape(chives) + tm_fill('trees_n', title = 'plant population',n=20, pal = "YlGn") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
legend.title.size = 0.9,
main.title = 'plant population, Chicago 2018',
main.title.size = 0.9)
print(plant_map)
#plant percentage map
plant_percentage_map = tm_shape(plant) + tm_fill('plant_percentage', title = 'percentage of plants',n=20,pal = "YlGn") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
legend.title.size = 0.9,
main.title = 'Percent population of plants, Chicago 2018',
main.title.size = 0.9)
print(plant_percentage_map)
#heat island index
heat_map = tm_shape(chives) + tm_fill('heatisl', title = 'heat island index',n=20, pal = "YlOrRd") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
legend.title.size = 0.9,
main.title = 'Heat island index, Chicago 2018',
main.title.size = 0.9)
print(heat_map)
#human map
human_map = tm_shape(chives) + tm_fill('acs_population', title = 'human population',n=20,pal = "RdPu") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
legend.title.size = 0.9,
main.title = 'human population, Chicago 2018',
main.title.size = 0.9)
print(human_map)
#human percentage map
human_percentage_map = tm_shape(human) + tm_fill('human_percentage', title = 'human percentage',n=20, pal = "RdPu") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
legend.title.size = 0.9,
main.title = 'Percentage of humans, Chicago 2018',
main.title.size = 0.9)
print(human_percentage_map)
#Spatial contiguity weight matrix PLANT
w.queen <- queen_weights(plant, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, plant["trees_n"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(plant), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}), border = "#333333", lwd=0.2)
title(main = "Plant population LISA (Queen)")
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")
lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(plant),
col=sapply(lisa_p, function(x){
if (x <= 0.05) return (p_colors[2])
else return(p_colors[1])
}),
border = "#333333", lwd=0.2)
title(main = "Plant Population LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")
#Spatial contiguity weight matrix HEAT
w.queen <- queen_weights(heat, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, heat["heatisl"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(heat), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}), border = "#333333", lwd=0.2)
title(main = "heat island index LISA (Queen)")
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")
lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(heat),
col=sapply(lisa_p, function(x){
if (x <= 0.05) return (p_colors[2])
else return(p_colors[1])
}),
border = "#333333", lwd=0.2)
title(main = "heat island index LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")
#Spatial contiguity weight matrix HEAT
w.queen <- queen_weights(human, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, human["acs_population"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(plant), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}), border = "#333333", lwd=0.2)
title(main = "pop density for humans LISA (Queen)")
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")
lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(heat),
col=sapply(lisa_p, function(x){
if (x <= 0.05) return (p_colors[2])
else return(p_colors[1])
}),
border = "#333333", lwd=0.2)
title(main = "pop density for humans LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")
#create dataset for plants and heat island index
hp = plant["plant_percentage"]
hp$heatisl = heat$heatisl
hp$pop = human$human_percentage
#correlation test with plant and heat island index
correlation = cor(hp$heatisl, hp$plant_percentage, method = "pearson")
head(correlation)
## [1] -0.08542622
#Overlay all the maps together
plant_bubble_map <- tm_shape(plant) +
tm_symbols(size = "plant_percentage",
col = "green",
alpha = 0.5)
heat_bubble_map <- tm_shape(heat) +
tm_symbols(size = "heatisl",
col = "red",
alpha = 0.5)
human_bubble_map <- tm_shape(human) +
tm_symbols(size = "human_percentage",
col = "blue",
alpha = 0.5)
all <- plant_bubble_map + heat_bubble_map + human_bubble_map
plantonly <- plant_bubble_map + heat_bubble_map
tmap_mode("plot")
## tmap mode set to plotting
print(plantonly)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
print(all)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
print(plant_bubble_map)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
print(heat_bubble_map)
print(human_bubble_map)
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.